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I. INTRODUCTION 

In the last decade enormous progress has been 
achieved in fabricating designed quantum systems on the 
nanoscale. Typical examples include the assembling of 
quantum dot arrays to form artificial molecules 0,0) the 
implementation of syntheziscd molecular structures in 
mesoscopic devices to develop molecular electronics 
and the tailoring of artificial atoms to obtain fully con- 
trollable quantum bits as the basis for quantum infor- 
mation processing |5(. A characteristic feature of such 
devices is that they consists of a relatively small number 
of relevant degrees of freedom, which are embedded in a 
much larger surrounding. 

Of particular importance are electronic transport pro- 
cesses, which are then determined by an intimate in- 
terplay between electron-electron correlations and in- 
teractions with the environment. In this context, the 
simplest model, well-known as the so-called spin-boson 
model @, Q, consists of a two-site chain where a sin- 
gle charge is transferred between the sites via tunnel- 
ing hybridization. The charge dynamics can be coherent 
or incoherent depending on the interaction with the en- 
vironmental degrees of freedom consisting e.g. of resid- 
ual vibronic and phonon degrees of freedom or electro- 
magnetic modes. In the past, the spin-boson model has 
been proven to capture essential properties of such elec- 
tron transfer systems including various sorts of quan- 
tum phase transitions. Together with its extensions to 
multi-site systems it describes a broad field of phenomena 
comprising e.g. quantum Brownian motion Q, Kondo 
physics (3 , Luttinger liquid behavior , atomic quantum 
dots |9j , and charge transfer in photosynthesis jlOj . How- 
ever, a fundamental limitation of the spin-boson model 
is its restriction to single charge transfer, which becomes 



particularly severe, when the number of charges and their 
dwell time on a certain site can be controlled. 

Recently, we studied a generalization to correlated 
charge transport based on a dissipative one-dimensional 
Hubbard model [l]]. In its isolated form (no dissipa- 
tion) the latter one has been extremely well analysed 
in various parameter domains and has served as a stan- 
dard model to reveal the nature of many-body corre- 
lations [H|- While those studies have been performed 
mostly in the thermodynamic limit to determine equi- 
librium properties and quantum phase transitions, the 
situation in tailored quantum systems is different: here 
the dynamics far from equilibrium and in presence of a 
dissipative surrounding is to be investigated, where typ- 
ically only few excess charges occupy a relatively small 
to moderate number of sites. One goal lies in the control 
of quantum properties in the transport, which inevitably 
happens to be influenced by energy exchange with and 
fluctuations from a phononic background. 

The nonequilibrium dynamics of correlated charge 
transport in dissipative environments is an extremely 
challenging task. This is mainly due to the fact that 
on the one hand away from the limit of very weak in- 
teraction with the phonon bath perturbative approaches 
cannot be applied and on the other hand at lower temper- 
atures phonon induced retardations become long-ranged. 
Hence, so far either completely coherent transfer (no or 
very weak dissipation) [Tj, [Tj| or completely sequen- 
tial transfer in the limit of very strong Coulomb repul- 
sion [lfj, [0| , which effectively reduces to a single charge 
problem, has been considered. 

In the last decade real-time path integral Monte Carlo 
techniques (PIMC) have been proven to provide an ef- 
ficient and numerically exact tool to simulate the time 
evolution of dissipative tight-binding systems in ranges 
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where perturbative schemes are prohibitive |7], Il7l | . The 
method has been substantially improved in the last years 
to tackle also large r sy stems (a larger number of sites) 
and many charges |18L Il9| . These simulations are chal- 
lenging in that the dimensionality of the Hilbert space 
to be sampled strongly increases with the number of 
particles, which in turn aggravates the so-called dy- 
namical sign problem 20]. The latter one originates 
from the oscillatory character of the time evolution op- 
erators generating "signals" due to interferences. We 
note in passing that imaginary time path integral Monte 
Carlo approaches have already been successfully applied 
to extract thermodynamic properties of the Hubbard 
model [U H3, HH- An alternative approach has been 
developed very recently by means of numerical renormal- 
ization group techniques (NRG) and its time dependent 
generalizations p3 |. which must be seen as being com- 
plementary to the PIMC scheme as the former ones are 
particularly adapted to the regime of very low temper- 
atures, while the latter one shows its best efficiency in 
intermediate temperature domains. 

The numerical results in Ref. elucidated the crucial 
relevance of symmetries of the full Hamiltonian compris- 
ing system, bath, and interaction. Namely, it has been 
shown that the generic dipole-dipole coupling between 
many-electron system and bath leads to the existence 
of a decoherence free subspace (DFS), i.e. a subspace of 
the total Hilbert space which is invariant under the time 
evolution. Interestingly, similar kinds of DFSs have been 
analysed in the context of quantum information process- 
ing p, however, in many-body systems their structure 
and impact on the transport dynamics is quite differ- 
ent. For instance, in contrast to the naive expectation 
for longer times the dissipative population dynamics of 
the many-body states approaches a steady state which is 
not of Boltzmann type and sensitively depends on the ini- 
tial preparation. In this paper we give a more elaborate 
account of this work and provide a complete analysis of 
the symmetry properties of dissipative one-dimensional 
Hubbard chains. In fact, a DFS turns out to be just a 
special kind of an invariant subspace embedded into a 
whole class of invariant subspaces. This characterization 
provides a detailed understanding of the relaxation to- 
wards a non-Boltzmann type of equilibrium and possibly 
opens new control mechanisms for fermionic transport. 
It further serves in the regime of completely incoherent 
transfer as a starting point to approximate the dynamics 
of the reduced density matrix by simple master equations 
with transitions rates obtained from golden rule calcula- 
tion, thus generalizing the rate concept known from single 
charge transfer ?j . Eventually, the existence of invariant 
subspaces may substantially improve the efficiency of the 
PIMC algorithm and particularly soothe the dynamical 
sign problem. 

The paper is organized as follows. We start in Sec. UTI 
with the description of the dissipative Hubbard model. 
In the next two sections, Sec. IIIII and Sec. IIVI the sym- 
metries of the model are analyzed in detail. This pro- 



vides in Sec. [V] insight into the dynamical and equilib- 
rium populations, before we give a brief description of 
the PIMC method and its description in terms of a rate 
model in Sec. lVIl Applications to specific realizations are 
discussed in Sec. IVIII At the end are short summary is 
given and some conclusions are drawn. 

II. DISSIPATIVE HUBBARD-MODEL 

The system is modeled by an open Hubbard chain with 
N sites of spacing qo 

N / II \ 

H S = £ U< CT ^ + y<AX-A- 

N-X 

+ E ^i{4,.d i+lta + h.c) , (1) 

where di yC7 and d\ are annihilation and creation opera- 
tors for electrons, respectively, with spin a on the site 
i. Ei are the bare energies of the levels, [/, denotes 
the strengths of the Coulomb interaction on site i, and 
A, are the tunneling matrix elements. The influence of 
the bosonic bath is given by a Caldeira-Leggett type of 
model |2jj, where the interaction with the bath 

is accomplished via a standard dipole coupling [13] 

Hi = -qoV £ c a X a + qlV 2 £ , (3) 

a a 

where 

N 

v= J2 [*-(^+i)/2]4x, ff ( 4 ) 

i— 1,<t 

is the polarization operator of the Hubbard chain. This 
way, the above model can be seen as a generalization of 
the spin boson model 0, of charge transfer to the case 
of many electrons and many sites. 

In the sequel we deal exclusively with distinguishable 
charges and thus, in the case of fermions, concentrate on 
the two particle/opposite spin sector, see Ref. ^l). We 
note in passing though that all results presented in this 
paper are fully transferable to any two particles as long 
as they can be physically distinguished. Based on this 
analysis, the generalizations to more charges/particles 
may be tedious in detail, but straightforward in prin- 
ciple. The case of indistinguishable particles is related 
to the bosonic Hubbard model and will be discussed in a 
subsequent paper. 

A suitable basis for the electronic degrees of freedom 
is now given in terms of the localized many-body states 
(LMBS) 

Sioc = {|s ; ,s t )a}_s< s i, s t<s , (5) 
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where \s^, s')a denotes an antisymmetrized state with 
the a = I (<t = t) fcrmion being localized on site i = 
s l +S+ 1 (i = s T + 5+1), where S = (d- l)/2, i.e. 

|aW>A = ^= (11,1,^)12,1,^) - 11,1,^)12,4,/)) , 

(6) 

with, e.g., |1, J,, /) denoting a state with the first particle 
being in the spin-down state, localized on site i = + 
5 + 1. Since only antisymmetrized states will be used 
throughout the rest of the paper, the subscript A will 
be omitted in the following. The time evolution of the 
LMBS populations then follows from 

P s i^(t) = (aVlpWkW), (?) 

where p(t) — tre{W(t)} denotes the reduced density ma- 
trix of the electronic system, obtained from the full den- 
sity matrix W(t) = exp{-iHt/H)W(0) exp(iHt/h) with 
H = Hs + Hb + Hi after integrating out the environ- 
mental degrees of freedom (see also Sec. IVI A"|) . 

As shown in 11], the reduced dynamics Q is strongly 
determined by the symmetries of the underlying Hamil- 
tonian of the total compound. In fact, the dipole type of 
system-bath coupling provides a symmetry which gives 
rise to a decoherence free subspace (DFS). The projector 
onto this subspace commutes with Hs and Hi, i.e. with 
the polarization V , meaning that it is not affected by the 
dissipative dynamics of the environment. Here we argue 
that such a DFS is actually just a special case of a more 
general property: The decomposition of the Hilbert space 
into invariant subspaces between which neither the free 
Hamiltonian Q nor the bath can induce transitions. Ac- 
cordingly, when populating initially only one of these sub- 
spaces, the dissipative system's dynamics remains com- 
pletely restricted to it. This has profound consequences: 
(i) Only subspaces of the full Hilbert space must be sam- 
pled in the Monte Carlo simulations which substantially 
reduces the dynamical sign problem and (ii) for appro- 
priate bath parameters the full quantum dynamics can 
be expressed into simple rate equations, which intrinsi- 
cally obey the corresponding symmetries and lead to the 
correct equilibrium populations. 

The analysis of the symmetry properties and its re- 
lation to invariant subspace is most conveniently done 
by exploiting that the rt-particle dynamics on a one- 
dimensional lattice can be mapped onto an effective sin- 
gle particle diffusion on a n dimensional one. For in- 
stance, as shown in Ref. [26(, for two distinguishable 
fermions on a one-dimensional lattice with d sites a 
mapping exists onto one particle dynamics onto a two- 
dimensional square lattice with d 2 sites, see Fig. ^ This 
mapping turns out to be a very convenient tool on 
the one hand to visualize the relation between many- 
body dynamics and dissipation and on the other it may 
serve as a starting point to apply pcrturbative methods 
as e.g. the non-interacting blip/cluster approximation 
(NIBA/NICA), which are known to be powerful means 
to capture the single particle motion. The rate equations 
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FIG. 1: 2d mapping of the LMBS © for d = 3. 

mentioned above are a direct consequence of this type of 
approximation. 

III. SPIN-PERMUTATION SYMMETRY 

Since the dissipative Hubbard Hamiltonian contains 
no terms depending on the fermions' spins, it is invariant 
with respect to any permutation of the latter. This leads 
to the most robust symmetry of the model, namely, 

[T, H] = [T, H s ] = [T, Hi] = [T, H B ]=0, (8) 

where we introduced the spin-permutation operator 

_F| S V> = |sV>, (9) 

which, for two particles, is equivalent to a simultaneous 
spin flip. 

Note that T commutes not just with the total Hamilto- 
nian H, but with Hq and Hi separately, such that we can 
easily find an electronic basis of simultaneous eigenstates 
of T and V . For example 

B mp ee B+ p UBfl ip ; 

B+ p ee {\aW = s L )} W < s 

L V2 J s i< s T 

Sflip = {7f(l si ' sT >-l sT ' si >)} si<sT ' ( 10 ) 

where B^ ip only includes eigenstates of T with eigenvalue 

±1, respectively. Transitions between and B^ ip are 
globally prohibited, independent of the system parame- 
ters d, Ei, Ui, A,-, or any bath parameters. Therefore, 
the electronic Hilbert space H decomposes into two dis- 
junctive invariant subspaces ee span{S^ ip }, i.e. 

H = g + ug-, g+ng- = 0, ng+ c g+ , hq- c g- , 

(ii) 
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with distinct routes to thermal equilibrium. As a conse- 
quence, the LMBS populations J7J) exhibit different relax- 
ation time scales and, even at high temperatures, a non- 
Boltzmann equilibrium distribution. Furthermore, the 
latter explicitly depends on the initial preparation. For 
example, for an energetically completely degenerate sys- 
tem (i.e. Ei — Ui = 0) the populations for the d(d+ l)/2 
states in B^ ip equilibrate towards 2Pq /d(d +1), while 
those for the d(d— l)/2 states in B^ ip equilibrate towards 

2Pq I ' d(d— 1), where Pq denote the initial populations of 
the subspaces G ± , respectively (Pq + Pq =1). Accord- 
ingly, even for sufficiently strong damping and/or high 
temperatures, the equilibrium populations of the local- 
ized states © do not coincide with a Boltzmann distri- 
bution. 

For the special case d = 2, it is straightforward to show 
that the basic spin-permutation symmetry fully explains 
the existence of a decoherence-free state as reported in 
Ref. ^l| (namely in form of the only element of B^ ip ). 
The robustness of this DFS with respect to changes in 
system and bath parameters directly follows from the 
robustness of the spin-permutation symmetry itself. For 
d > 2 wc postpone a more detailed discussion and first 
analyse in the next section additional symmetries related 
to particular parameter ranges for Ei,Ui, A,. We already 
note, though, that the spin-permutation symmetry can 
be destroyed e.g. by applying external magnetic fields, 
which in turn would also allow to control its influence on 
the many-body dynamics. 



IV. PSEUDO-ANGULAR MOMENTUM BASIS 
AND INVARIANT SUBSPACES 

For the tight binding lattice with the localized ba- 
sis 101, it is very suggestive to identify the discrete posi- 
tions of the charges and with the discrete eigenval- 
ues of the z-components JJ and Jj , respectively, of two 
pseudo-angular momentum operators. Accordingly, the 
latter ones obey J 12 = J T2 = h 2 S(S + l) with d = 25 + 1. 
For single charge transfer on two sites this leads to the 
famous spin-boson model where the polarization is then 
proportional to the Pauli matrix a z . Now, for two charges 
the polarization V is proportional to + Jj . Hence, for 
the analysis of symmetry properties it is convenient to 
introduce an alternative basis of the fermionic Hilbert 
space following from a total pseudo-angular momentum 
J = J* + JT . The consequences of this formulation are 
studied in this section. 



A. Collective states 

Within the pseudo-angular momentum representation 
the localized states in B\ oc (JSJ) are simultaneous eigen- 
statcs of J^ 2 , Jj-, J T2 , and J\ due to 

|s^, s^) = \j^ = S, j 1 ' = S; m)- = s^, toj = s') , (12) 



while, as already mentioned, the polarisation operator 
can be expressed as 



V = 



Jl + J I 



h 



(13) 



Here J z denotes the z component of J = J* + . An 
alternative electronic basis set is then given by the simul- 
taneous eigenstates of these two operators, i.e., 



S.7 = {IV'm )}0<7<d-l,-i<m<j 



(14) 



where J \ if). 

l(J) 



H 2 j{j + and J z \^>) 



(i)\ 



Km\ipm )■ Note that, like the localized states in B\ oc , the 
collective states in Bj are also eigenstates of the polariza- 
tion operator V (|13[) . as well as of the spin-permutation 

operator T © (^|Vm) - |V>m))- 

The transformation between the two basis sets is con- 
veniently given in terms of the (real-valued) Clebsch- 

Gordan coefficients H3 (s^s^m) = (5, S, s l , s T |^ } ). 
Hence, defining 

C^^VIV^t) (15) 
and exploiting (s^, \tjjm) = for m ^= s^+s^ , we obtain 



m.in{S,m-\-S} 

y 

Z ^ s^,m — s 

s^— max{ — S,m — S} 

d-1 

= E C V )S tlV' s i +s T) . 
j=| s l+ s T| 



s+ .m — st 1 ' ' J 



(16) 



With the aid of Eq. ED on e can now easily express the 
LMBS populations (Q) in terms of the corresponding ones 
in the pseudo-angular momentum basis, 



yielding 



d-1 . 2 

j=| s J-+s T 



(17) 



d-l 

£ 

j,j' = \s l +s r \ 



C« sT Cif] T Re{(^'| sT |K*)l^l s t)} 



(18) 

Equation (|18|l also governs the transformation of the ini- 
tial populations between the two basis sets. Particularly 
simple expressions are gained for an initially factoriz- 
ing system-bath state where, as in many experiments, 
the system is prepared to reside on some localized state 
|sg,sj). The corresponding initial populations -P,-, m (0) 
are then obtained from 



Pj.m (0) 



Cj't if m=sJ + 4, j > 



else . 



(19) 
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B. Invariant subspaces 

The formulation based on the total pseudo-angular mo- 
mentum now allows to identify a further symmetry of the 
dissipative Hubbard model. For this purpose we first con- 
sider the case, where the full Hamiltonian H commutes 
with the total pseudo-angular momentum J 2 , thus lead- 
ing to invariant subspaces similar to the symmetry noted 
in Sec. IIIII Classically speaking, H preserves the angle 
between the individual pseudo-angular momenta and 
JT. 

Accordingly, we look for a representation of the elec- 
tronic Hamiltonian Hg in terms of the components J X1 
J y , J z . Turning first to the nearest-neighbor hopping 

term Hg = J2i,a Ai(<4 CT d l+l rT + h.c), it is straightfor- 
ward to show that with respect to B\ OCl it exhibits the 
same structure as J Xl namely, 

(^,sT| jff A| s l )S T) = ^ (^,|T|JzkV) =0, 

provided all couplings are non-zero. 

Since the non-vanishing matrix elements of J x are 
given by 

± 1, s T | J,|s V) = yS(S + 1) - si(si ± 1) (21) 

(and (s^,s^ ± 1| J x |s , s*) accordingly), one finds that 
H{g oc J x is always true for d = 2. For d > 2 this holds 
only if the tunnel couplings (hopping amplitudes) are 
chosen properly, e.g. for d = 3 one must have Ai = A2, 
i.e. isotropic coupling, and d > 4 requires quite specific 
values for the coupling strengths. We will discuss this 
point in more detail below. 

Next, we turn to the on-site terms in Hs comprising 
the on-site energies and the Coulomb interactions. Since 
the latter one effectively gives a contribution to the for- 
mer ones, we find that if Ei and Ui, i = 1,...N, are 
distributed up to an overall factor according to the ma- 
trix elements of J z , i.e., 

(s^V.kV) = (s l + s t )<Sjw<Sjt= s t , (22) 

then the electronic Hamiltonoperator can be cast into the 
form 

H s = EJ Z + AJ X , (23) 

with appropriate constants E and A. It then follows 
immediately that 

[H, J 2 ] = [H S ,J 2 } = [Hi , J 2 ] = [Hb, J 2 ] = . (24) 

The total pseudo-angular momentum is thus a con- 
served quantity with respect to the full dissipative sys- 
tem. Equivalently, for a dissipative Hubbard model with 
Hg of the form (12;j[) both the free system as well as 
the bath can induce only transitions between states with 
the same total pseudo-angular momentum. Of course, 



Eq. (|23(l is a sufficient but not necessary condition for 
Eq. (|24|l to hold. Obviously, any Hs and Hj where the 
system part can be written in terms of J x , J y , and J z 
(and powers of these operators) fulfill Eq. I|24|l : in partic- 
ular, the coupling between system and bath may also be 
nonlinear in the system operator. Indeed, below we will 
present an example where Eq. (12.'j[) does not hold, while 
the latter relation still applies. 

With respect to the reduced electronic dynamics, 
Eq. (P4~|l leads to consequences similar to those described 
in Sec. IIIII First, the electronic subspace TL of the full 
Hilbert space can be further decomposed into disjunctive 
invariant subspaces, 

n = u^tcfe , (25) 

where Qj denotes the subspace spanned by all eigenstates 
of J 2 with a fixed total pseudo-spin j, 

^span^j),...,^)}, (26) 

with the properties 

Q k DGj=0 for all k j= j , HQj C Q j . (27) 

To put it differently, with the above decomposition we 
have identified the irreducible representations of the sum 
of two angular momentum operators each with angular 
momentum S = (d — l)/2. Unlike the spin-permutation 
symmetry, now the number of invariant subspaces equals 
d, thus explicitly depending on the length of the tight 
binding lattice. Furthermore, the structure of these sub- 
spaces turns out to resemble linear chains since due to 

(^£\H S \^) = ES m , =m 

+ y (VJ(J + 1) - Mm + l)<5 m / =m+ i 

+ Vj(j + 1) - m(m - l)<5 m '= m _i) • (28) 

Transitions thus are only possible between nearest neigh- 
bors with respect to m between states \ipm^) with fixed 
j- 

To summarize the main outcome of this analysis, we 
have shown that an electronic Hamiltonian obeying ill'-'il) 
leads to a Hilbert space basis which can be decomposed 
into one singlet part, one triplet part, one quintet part, 
. . . , up to a 2(d — 1) + 1 part, each of which has a lin- 
ear structure and between which transitions induced by 
the Hamiltonian are not possible whatsoever (see Fig. [21 . 
The subspaces solely interact via the coupling to a com- 
mon bath. The singlet space is exceptional in that it is an 
eigenstate of Hj with vanishing eigenvalue for the polar- 
ization operator, thus being protected by any dynamical 
time evolution of the dissipative system, i.e. a DFS. 

V. POPULATIONS OF THE LOCALIZED 
MANY-BODY STATES 

For transport processes the populations of the LMBS 
are of particular relevance. Here, we show how the exis- 
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FIG. 2: 2d mapping for d — 3 and the angular-momentum 
states 1141 . Black lines indicate the couplings for ei = U — 0; 
the thickness of the lines is proportional to the square of the 
coupling strengths. Gray solid lines indicate additional cou- 
plings for U 7^ 0, tj 7^ 0; gray dashed lines indicate additional 
couplings for anisotrope couplings Ai 7^ A2. 



tence of the ideal symmetry l|24fl may simplify their cal- 
culation considerably. This gives also direct access to the 
stationary equilibrium populations which are approached 
for sufficiently long times. 



A. Symmetries and population dynamics 

Due to the disjunct structure of the full Hilbert space, 
see Eq. (|2T|) . not only the total population is conserved, 
but also the individual ones within each invariant sub- 
space. In particular, for a factorizing initial preparation 
W(0) = po exp(— /3Hb), where po = \(po){\<po\ presents a 
projector on some pure fermionic state \ipo), subspaces 
with no overlap with |<po) do not participate in the dy- 
namics at all. This leads to a reduction of the dimension- 
ality of the relevant Hilbert space and thus of the com- 
plexity of the system to be treated numerically. If the 
initial electronic population is even restricted to a single 
subspace, i.e. \ipo) £ Gj, the corresponding dynamics is 
identical to that of a single particle moving on a linear 
chain with just 2j + 1 sites. For \ip ) £ Bi oc , this sce- 
nario only applies to the states \—S, —S) (i.e. |^_/^_]\)) 

and \S,S) (i.e. \ip d d ~^)) with d = j + 1. We note in 
passing that this property also applies to the situation 
|yo) £ Gj U Go since Go does not exhibit any dynami- 
cal evolution at all. For other initial preparations within 
B\ oc , however, direct phonon induced couplings between 
the participating subspaces mediated by the common en- 
vironment 26] have to be taken into account, which do 
influence the dynamics but conserve the population in 
each subspace. For initial states which have overlaps 
with all invariant subspaces (i.e. \—S, S) or IS*, —S}), the 
dynamics only takes places in a (d 2 — 1) dimensional sub- 
space since again Go does not participate. 



Due to the invariance of the subspaces Gj with respect 

to H, the cross terms (^.Lf IpWIV'H+st) in Ec L- UHl are 
only non- vanishing if the initial state \ipo) has an overlap 
with both Gj and Gj' ■ For most preparations in a LMBS, 
this significantly simplifies Eq. I|18l) . This holds especially 
for preparations where both fermions enter at the same 
end of the chain (i.e. |— S, —S) and \S, S), both £ Gd-i), 
when all localized populations can be expressed in terms 
of the pseudo-angular momentum populations alone. For 
all other localized preparations, however, the cross terms 
play a crucial role in the short-to-intermediate time do- 
main. On the other hand, using Eq. (|18l) to obtain in- 
formation about these cross terms also eludicates the dy- 
namical evolution of coherences in the system. 



B. Equilibrium populations 

For long times a dissipative system reaches a stationary 
state, typically a thermal equilibrium where all states are 
occupied according to a Boltzmann distribution. In the 
case studied here, namely, a dissipative Hubbard model 
with the symmetry relations Ij24(l . the same applies not to 
the total Hilbert space, but rather its disjunct subspaces. 
In particular, for a state \iprn) € Gj one obtains 



P°° m = lim P jtm (t) = P c 



-H/3e { ri 



y 3 > 



(29) 



where the initial population in Gj is P 
Ei=-i4m(0), and hett = {^\H S \^). Accord- 
ingly, the final occupations of the states \tprn) do not 
obey a Boltzmann distribution, i.e., 



P 



1 1 \ 



poo 
j',m' 



(30) 



The same is true for the states in B\ oc . Since in thermal 
equilibrium all cross terms in Eq. I |18| l vanish (P°? sl — 



poo 



we have 



d-i 



(j) 2 p o 



poo _ \ " Mj) 

j=\sl+sT 



U) 



si + s 



(31) 



with Zj = £ J m , = 



. Apparently, the equilibrium 



populations strongly depend on the initial preparation, 



i.e. on P, 







VI. 



QUANTUM MONTE-CARLO 
SIMULATIONS 



So far we have focused on a convenient representation 
of the states within the fermionic subspace of the to- 
tal Hilbert space. To obtain the populations Eqs. Q 



7 



and l|17|) explicitly, however, one must find an adequat 
treatment of the environmental degrees of freedom as 
well. The path integral approach has been proven to al- 
low for an exact elimination of the latter and to provide in 
combination with Monte Carlo techniques a numerically 
exact evaluation. 



A. Path integral representation 



The expectation value of a system observable A 



(A(t)) = T rs {p(t)A} 



(32) 



requires the knowledge of the reduced density operator 

p(t) = Tr B {e- lHt/h W{0)e lHt/h } . (33) 

Here, H denotes the full Hamiltonian 0-©, W(0) is 
an initial density, and the trace is performed over the 
bath degrees of freedom only. The details of the path 
integral representation for p(t) have been given else- 
where Here, the standard expression known 
e.g. for the spin-boson model has to be extended to cap- 
ture the dynamics of two distinguishable fermions. Ac- 
cordingly, the reduced density operator is expressed as a 
double path integral along a Keldysh contour with for- 
ward s 17 and backward s" paths corresponding to the ba- 
sis set B\ oc . The impact of the dissipative environment 
appears as an influence functional introducing arbitrar- 
ily long-ranged interactions in time along and between 
the paths. It is convenient to switch to the sum and 
difference coordinates r) a ~ s a + s CT and £ CT ~ s a ~ s CT , 
respectively, so that one arrives for the expectation value 
at the exact expression 

(A(t)) = Ivff lvZa%&K%& e X p{-$[77,|*]} , 

(34) 

where a[ff, £] is the measurement functional correspond- 
ing to A in terms of the combined system paths fj{t) — 
(rf(t),r] l (t)) and^) = (^(t),^(t)). Furthermore, K, is 
the bare action factor in absence of a reservoir, and the 
influence functional reads 



ds du ■ 
o Jo 



{[|*(s)-e] L'(s-u) [i(u)-e 



+ i 



f(s) • e L"(s - u) 



£(s) • e rf(s) ■ e 



%u) • e | 

(35) 



withe = (1, 1). The kernel L{t) = L'(t)+iL"{t) is related 
to the force-force auto-correlation function of the bath 
and is completely determined by the spectral density of 
its modes HA7jl. Further, \i = lim^^o h/3L(0). Note that 
even in absence of Coulomb interaction the two charges 
are correlated due to the coupling to the common heat 
bath. 



An analytical treatment of the expression l|34|l is in 
general not feasible, mainly due to the retardations in the 
influence functional which grow with decreasing temper- 
ature, roughly as Hf5. In this situation PIMC methods 
have been shown to be very powerful and numerically 
exact means to explore the non-perturbative range. Ap- 
pendix ^ explains in some detail our approach which is 
straightforwardly obtained from the formally very simi- 
lar case of a single dissipative particle ja, Q, [111 uM an d 
which can be easily extended to the cases of more than 
two or indistinguishable particles p6l |. 

A striking difference to the single-particle case, how- 
ever, is the existence of symmetries as described above 
and the consequent decomposition of the full system's 
Hilbert state into invariant subspaces. To exploit these 
symmetries for our numerical studies, we express Eq. I|34|) 
in terms of system paths J(t) and Mj(t), referring to the 
states IV'm^) <|14|) rather than in terms of s^(t) and (t). 
It turns out that the environmental influence can again be 
summarized in an influence functional which is obtained 
from Eq. 135fl by simply replacing ff(s) • e and £(s) • e with 
7] M = Mj(t) + Mj{t) and £m = Mj(t) - Mj(t), respec- 
tively, where Mj(t) and Mj(t) again denote paths on the 
forward and backward part of the Keldysh contour, re- 
spectively. The corresponding influence functional then 
reads 



$[vm,£m] 



t , s 

ds du[£ M (s)L'(s - u)£ M (u) 
o Jo 

+ i£ M (s)L"(s - u)t]m{u)} 

+ i- / ds£, M (s)r) M (s) ■ 
1 Jo 



(36) 



Note that for fixed J this is just the influence functional 
of a single dissipative particle residing on 2 J + 1 discrete 
states. 

Now, if Eq. 12411 holds, J{t) is conserved upon prop- 
agation with exp(±itH/h) (see Eq. IjAHl ). A change in 
J(t) can thus only be mediated by the initial prepara- 
tion Tr{iy(0)} and the measurement a[rjM,S,M\ at time 
i; otherwise J(t) remains constant. In terms of paths of 
J and M, Eq. (|34l) therefore eventually becomes 

(A(t)) = A 'ujM, 

Ji,.h=0 

x ^JuJ 2 [vm,£.m] exp{-$[?7M, 6/]} ,(37) 

where J\ and Ji denote the values of the piecewise con- 
stant J paths along the forward and backward parts of 
the Keldysh contour, respectively, and aj t ,j 2 [?7m, £m] is 
the measurement functional of the operator Pj^P^, 
with the projector Vj onto the subspace Qj. 

Further simplification can be achieved by exploiting 
that Aj lt j 2 (t) = if any of the subspaces Tre{ W(0)}Gj 1 , 
TibIW (0)}G ,j 2 or AQj 17 AQj 2 are empty. For example, 
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if the system is initially prepared in the localized state 
|— S, — S), one simply obtains 

(j4(t)>=4j-i,d-i(t). (38) 

Accordingly, in cases where the initial preparation does 
not impose transitions between subspaces with different 
J, the dissipative two-particle system can be decomposed 
into d independent dissipative single-particle systems. For 
PIMC simulations this presents a major progress since 
each of the systems representing the Aj lt j 2 (t) contribu- 
tions exhibits a significantly lower dimensionality than 
the full d 2 dimensional two-particle system. Hence, by 
evaluating Eq. I|37jl rather than Eq. I|34|) . the dynami- 
cal sign problem |20j | , which stems from interferences be- 
tween different quantum paths and grows exponentially 
with the size of the system, can be greatly reduced. In 
the quantum regime this allows for substantial compu- 
tational savings, despite the fact that the evaluation of 
Eq. I|37[) requires d(d — l)/2 independent PIMC simula- 
tions (note Aj 2jl (t) = A* Ju j 2 {t)). 

In case that Eq. 1|24[1 does not apply, an expression 
similar to Eq. (|37l) can still be obtained exploiting the 
spin-permutation symmetry. Therefore, the arguments 
presented above concerning the computational costs of 
PIMC simulations still apply, albeit to a lower degree. 

For further details of the approach and particularly of 
the strategy to soothe the d yna mical sign problem we 
refer to the literature 0, Il8l.ll9j . The basic ingredients 
are this: (i) One exploits the fact that the influence func- 
tional depends only linearly on the quasi-classical paths 
fj{t), so that in Eq. i|34|) the corresponding summations 
can be expressed as a series of simple matrix multiplica- 
tions and therefore be carried out explicitly; (ii) for the 
sampling procedure one chooses an MC weight with no 
long-time retardations. Since they are fully taken into 
account in the final accumulation process evaluating the 
exact expression (|34|) . the numerical exactness of the MC 
scheme is not impaired. This way, one achieves a strong 
decoupling of quasi-classical (if) and quantum (£) coordi- 
nates, which allows to store products of short time prop- 
agators independent of the MC-sampling. Eventually a 
speed-up with respect to the original method by a 
factor of about 100 is gained. Of course, this strategy 
can be applied for evaluating both Eq. (|3*4"|) as well as 
Eq. P7| l. 

B. Rate model 

For one-dimensional tight-binding systems a descrip- 
tion of the single particle population dynamics in terms 
of simple master equations with transition rates gained 
from golden rule calculations is known to be quite accu- 
rate for sufficiently strong dissipation and/or high tem- 
peratures (incoherent dynamics) |7lll9l|. Then, the exact 
dynamics approximately obeys 

^ = HP(t), (39) 



where P collects the single particle populations and the 
matrix R contains the transition rates between adjacent 
sites. These rates obey detailed balance reflecting the ex- 
istence of a thermodynamic equilibrium approached for 
long times, where the populations are Boltzmann dis- 
tributed such that RP(i — > oo) = 0. However, as we 
have shown in the previous sections, for the many-body 
time evolution this is no longer true, and formulating 

Eq. with P(t) = (P_ s ,-s(i),P-S,-S+l (*),-••) and 

the corresponding golden rule rates in R must inevitably 
fail to reproduce the correct dynamics. 

The idea is thus to start from the pseudo-angular mo- 
mentum basis (|14f) and to employ separate rate descrip- 
tions for each invariant subspace with initial and equi- 
librium populations according to Eqs. (|19|l and l|29|) . re- 
spectively, and with transition rates chosen according to 
known expressions like golden rule formulae 

MM The 

LMBS populations can then be obtained from Eq. l|18fl . 
This approach, however, only applies for LMBS pop- 
ulations for which the cross terms in Eq. Ijl8|l vanish, 
i.e. if the initial density matrix p(0) includes no finite off- 
diagonal elements with respect to the basis Bj. In terms 
of the path-integral picture, this restricts a rate approach 
to populations whose path-integral expression (|3TJl col- 
lapses to a single A,j 1 _j 1 (t). For instance, for a sys- 
tem obeying Eq. I|24|) and being prepared according to 
P-S,-s(0) = 1 or P S ,s(0) = 1, this holds for all LMBS 
populations, while for Po,o(0) = 1 only some of them 
can be reproduced this way (cf. Figs. |3]and0J). Note 
that Eq. 124|) is no necessary condition for this approach, 
as we could confirm for a system with d = 3, ti = 0, 
Ai = Aa, U = 5Ai, and P_x x(0) = 1, where despite a 
coupling between Go and Qi all LMBS populations with 
s l _|_ s 1 could be reproduced (not shown). 

If, tough, off-diagonal terms are present in p(0), they 
by no means can be captured by a simple rate formalism, 
and the dynamics of the corresponding LMBS popula- 
tions will, at least for short to intermediate times, escape 
a description along the lines of Eq. H39f) . This reflects 
the fact that, although unable to introduce transitions 
between them, the phonon-induced coupling still corre- 
lates the dynamics on otherwise disjunctive subspaces, 
an effect clearly beyond the scope of rate models. To 
what extent the presence of non-diagonal elements of the 
density matrix do influence the dynamics in specific sit- 
uations, however, necessitates a deeper analysis based on 
the exact path integral expression (|34|l and a generaliza- 
tion of what is known about standard spin-boson type 
of models 0. Nevertheless, we heuristically found that 
in the cases we investigated, LMBS populations which 
cannot be solely expressed in terms of pseudo-angular 
momentum populations still can be well approximated 
by 

P.i,.r(t) = Psi, s i (t) + (1 - e" 7 *) (P s T, sT " i*5>) , 

(40) 

where P s i s t(t) denotes the corresponding LMBS pop- 
ulation as obtained from Eq. (|39|l . and 7 is the small- 
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est finite modulus of the eigenvalues of the rate matri- 
ces of all invariant subspaces Gj which participate in the 
dynamics (cf. Fig. 0J. Whether Eq. H4()Jl also covers a 
more general set of parameters, however, lays beyond 
our knowledge as well as an explanation for observation 
that in cases where no finite non-diagonal elements are 
present, Eq. I|4U|) performs rather poorly. Thus, Monte 
Carlo simulations are of particular importance to obtain 
insight into the nonequilibrium dynamics of correlated 
charge transfer in cases where symmetries are broken. 



VII. APPLICATIONS 

In this section we apply the general results of the pre- 
vious sections to specific models, particularly, to tight- 
binding systems with d — 3 and d — 4 sites. The sim- 
plest case d = 2 is not explicitly addressed since it has 
been extensively studied in Ref. [llj . We only mention 
again that for this generalization of the ordinary spin- 
boson model, a decomposition in Q + = Gi and Q~ = Go 
exists for arbitrary system parameters, where the singlet 
state in the latter subspace constitutes a DFS. 



A. The case d = 3 

As already stated above, for isotropic couplings and 
all on-site energies and the Coulomb energy set to zero, 
one has a full decomposition of the Hilbert space. Cor- 
responding data for the dynamics of the localized pop- 
ulations (0) are depicted in Figs. and Q] for an ohmic 
spectral density with exponential cutoff, 

J(lo) = 2naLuc- uj/uJc , (41) 

where a = 0.25, uj c = 5A, and h(3 = O.lfiA (the bath 
setup is the same for all populations shown here). Note 
the change in the equilibrium populations due to the dif- 
ferent initial preparations. Turning on the Coulomb 
interaction, i.e. U ^ 0, a coupling between {ip^} and 
|4 2) ) is introduced (cf. Fig. EJ), with (^ 2) |#s|V>q 0) ) = 
(\/2/3)[/, such that only two invariant subspaces remain, 
i.e. 

U d=z = £ou2 U Gi , (42) 

with C7ou2 = <?oU<?2- However, when fixing the interaction 
strength to the specific value Uc = £-1 + £i — 2e , a full 
decomposition is recovered, see Fig. |SJ 

Note that, without an external magnetic field, even for 
arbitrary system parameters the particle-exchange sym- 
metry cannot be broken: Albeit \ipQ } then typically cou- 

(21 (21 (2) 

pies to \ip_'l), \tp ), and \tp\ ), no couplings between 
angular-momentum states with even and odd J exist, 
such that always a non-Boltzmann distribution for the 
equilibrium populations is obtained (cf. Fig.EJ). 




A t 

FIG. 3: Initial dynamics of a two-particle system for d = 
3 with completely degenerate onsite energies and coupling 
strengths Ai = A 2 = A and Uc = 0. Shown are PIMC data 
for the localized populations Eq. J7J as obtained from Eq. I34H 
(empty black symbols) and Eq. I37H (filled gray symbols); 
gray lines denote results from the rate approach described in 
Sec. I VI 51 Shown are P_i,_i (circles), P-i,o (triangles up), 
f— i,i (triangles left), Po.-i (triangles down), Po,o (squares), 
Po.i, ("+"), Pi-i (triangles right), Pi, ("x"), and P M (di- 
amonds), with equilibrium values of 1/5, 1/10, 1/30, 1/10, 
2/15, 1/10, 1/30, 1/10, and 1/5, respectively. The inset ex- 
hibits are more detailed view of the populations with no initial 
population. 




A t 

FIG. 4: Same as Fig. 01 but with initial preparation on |0, 0). 
Shown are P-\,-i (triangles left), P_i,o (squares), P_i,i (tri- 
angles up), P ,-i ("+"). Po.o (circles), P ,i, ("x"), Pi _i 
(triangles down), Pi ; o (diamonds), and Pi.i (triangles right), 
with equilibrium values of 2/15, 1/15, 2/15, 1/15, 1/5, 1/15, 
2/15, 1/15, and 2/15, respectively. Dashed gray lines denote 
P_i,i(t), Po,o(t), and Pi,-i(t) EJ. 



B. The case d = 4: Example for a partial 
decomposition 

For a tight-binding lattice with d = 4 sites, a full de- 
composition of the Hilbert space is already not possible 
for vanishing energies and isotropic hopping terms. Nev- 
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FIG. 6: Same as Fig. GO but with A 2 = 0.7Ai, ei = -1.65Ai, 
£2 = 0, £3 = 3.5Ai, and 17c = 5.1, with equilibrium values of 
0.1884, 0.1330, 0.0812, 0.1330, 0.1567, 0.0795, 0.0812, 0.0795, 
and 0.0673. Gray lines here denote equilibrium populations 
as obtained from a Boltzmann distribution if all symmetries 
were broken. 



ertheless, one can still find a partial decomposition 



7~ld= 



Go U g 2 n G 



1U3 • 



(43) 



where C7 1U3 = Si U S3, and So, S2, and Gim are again 
disjunctive and invariant with respect to H (see Fig. 01 . 
A Boltzmann distribution again can be observed for the 
equilibrium distribution of \ipm ) states for each invariant 
subspace but not for the total Hilbert space. Contrary 
to Si and G3, a combination of invariant subspaces like 
^iu3 w ill hi general not exhibit a linear structure with 
respect to the tunnel couplings. For a energetically non- 
degenerate systems, further couplings between So and G2 
are introduced. However, even for arbitrary on-site en- 
ergies Ei, Coulomb energies U, and tunnel couplings Aj 
the decomposition of the Hilbert space can never be com- 
pletely destroyed due to the underlying spin-permutation 
symmetry discussed in Sec. IIIII 



C. Quasi-invariant subspaces 

For systems with a high symmetry (e.g. when Eq. I)24|) 
applies), even a minute change in the system's param- 
eters can lead to a symmetry breaking. Along goes a 
modification of the structure of the Hilbert space: For- 
merly invariant subspaces merge such that their number 
decreases. Since the size of the new invariant subspaces 
exceeds that of the former ones, the equilibration pro- 
cess within such a larger subspace is expected to slow 
down accordingly. If, however, the coupling between 
formerly invariant subspaces remains rather small, the 
relaxation process may in fact occur on two (or more) 
time scales: One (or more) governing the internal dy- 
namics within each former subspace and another one, 
strongly separated from those, governing the overall ap- 
proach to equilibrium. In this situation, one may still 
speak of quasi-invariant subspaces which establish quasi- 
stationary states on an intermediate time scale. 

For example, a weak particle interaction with Uc — 
0.05A in an otherwise energetically degenerate system 
with d = 3 leads to a weak coupling between IV'o^) and 
|V>o ) (cf. Fig.[5J. Consequently, instead of the formerly 
invariant subspaces Go, Si, and S2, now only Si and 
Gou2 = Go U G2 remain invariant with respect to H . The 
smallest eigenvalue of the according rate matrix R (see 
Sec. IVIB)l . describing the long-time dynamics [l9|. then 
becomes 0.010A, suggesting an equilibration timescale of 
the order of 100A -1 . The second smallest eigenvalue, 
however, exceeds the latter by a factor of more than 
20. Accordingly, for 1 < At < 100, the system will be 
trapped in a 'quasi-equilibrium', which still reflects the 
full symmetry and therefore significantly differs from the 
true equilibrium state (cf. Fig.[SJ. From an experimental 
point of view, this quasi-equilibrium might easily be mis- 
taken as the true equilibrium since on the one hand, the 
equilibration timescale can take very large values for sys- 
tems with only slightly broken symmetries, while on the 
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FIG. 8: Same as in Fig. El but with Uc = 0.05A. The gray 
dot-dashed lines show the equilibrium populations Pfj^j = 
Pi" = 0.1663 (circles and diamonds), P!°i i0 = P °°-i = 
P(U = Pifi = 0.0835 (triangles up, triangles down, "+", and 
"x"), 0.0834 (triangles left and triangles right), and 0.1666 
(squares), respectively. Solid dark-gray dotted lines represent 
the equilibrium populations for a corresponding system with 
U c = 0. 



other such a symmetry breaking is rather hard to iden- 
tify since it can easily occur at rather 'harmless' system 
parameters. Note that this phenomenon can also occur 
when applying a weak magnetic field which breaks the 
spin-permutation symmetry. 



D. Particle interaction 

An unique feature of many-particle systems is the in- 
teraction between its constituents, where the Coulomb 
interaction is the most prominent example. In our 
model CQ~©i Coulomb interaction between particles 
manifests itself as change of the energies of the local- 
ized states and influences the dissipative system in quite 
different ways: On the one hand, it can facilitate or break 
the symmetries discussed in Sees . Hill and II V| and accord- 
ingly the decomposition of the full Hilbert space into 
invariant subspaces, as has been intensively discussed 
above. On the other, as has been already pointed out by 
R. Marcus [2^, a different energy offset between adjacent 
sites A.Ej = Ej — Ej + \ introduces a phonon activation 
barrier between these sites, which is a unique signature of 
the interaction between system and phonon bath. It re- 
flects the fact that for a transition between adjacent sites 
energy fluctuations are necessary for the reorganization 
of the bath degrees of freedom after the tunneling of an 
electron. This latter process can exhibit a strong impact, 
since for larger Coulomb interaction Uj, the phonon ac- 
tivation barrier exceeds the electronic offset AEj by far. 
In fact, it can lead to rather counter-intuitive dynam- 
ics: If, for example, two electrons are initially prepared 
on the same site of the chain, any Coulomb interaction 
is expected to speed up the initial state's depopulation. 




Af 

FIG. 9: APo,o(t) for a system with d — 3, 6; = 0, and 
Uc I A = (circles), 2.5 (squares), 5 (diamonds), 7.5 (triangles 
up), 10 (triangles down), 15 (triangles left), and 20 (triangles 
right), initially prepared in |0,0). 

However, in competition to that the phonon interaction 
also generates an activation barrier, which in turn slows 
down any transitions between the initial and neighbor- 
ing states. Therefore, only sufficiently weak Coulomb 
repulsions will indeed lead to a faster emptying of the 
initial state; large repulsion strengths, in contrast, cause 
a slow-down of the initial depopulation dynamics and 
can even lead to a stabilization of the initial state. Fig- 
ure El visualizes this by comparing the initial dynamics 
of systems with varying interaction strengths Uc- Note 
that the turnaround takes places at around Uc — 2.5/iA, 
which coincides with the classical reorganization energy 
hh.W = (1/tt) J^duj J{lo)/lo. 

We note in passing that a quantitative discussion of the 
effect of phonon-induced interactions via the common en- 
vironment is, unlike to the case of two indistinguishable 
particles, not easy: While in the latter case one could 
simply compare the dynamics of two particles on a Id 
chain with Coulomb repulsion Uc with the QMC data 
for the dynamics of a single particle on a corresponding 
2d lattice 26], for two distinguishable particles this ap- 
proach is of no avail due to the intrinsically completely 
different structures of the Hilbert spaces. 

VIII. CONCLUSIONS 

In this paper we have analysed the correlated nonequi- 
librium dynamics of two interacting distinguishable 
fermions on a one-dimensional tight binding lattice em- 
bedded into a phononic environment. As a major result, 
we have revealed the crucial role of symmetries and pro- 
vided their full characterization. In particular, there is: 
(i) a spin-permutation symmetry which is fundamental 
in that it is conserved independent of system and bath 
parameters; it leads to a separation of the full Hilbert 
space into two disjunct subspaces. (ii) A conservation 
of a total quasi-angular momentum J corresponding to 
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[H, J ] = for certain topologies of system parameters. 
This symmetry determines a set of J = d (d number of 
sites) disjunct subspaces. The dissipative many-body dy- 
namics happens to occur in each of these subspaces inde- 
pendently with transitions between them only mediated 
by the initial preparation and/or the measurement ob- 
servable. Recently, similar kinds of symmetries have been 
discussed for quantum state transfer along networks as 
realized e.g. in linear arrays of Josephson junctions |29| . 
In any case, the time evolution in presence of the heat 
bath is not ergodic moving towards a stationary state 
which is not Boltzmann like distributed. Of particular 
relevance is the DFS, a subspace with no coupling to 
the bath at all, the initial state of which is preserved 
for all times. For d = 2 the nature of the many body 
DFS is equivalent to the DFS discussed in the context of 
quantum information processing for interacting spin-1/2 
systems, but differs from the latter for all d > 2. 

A direct practical consequence of the symmetry prop- 
erties is that the dynamical sign problem notoriously oc- 
curring in real-time quantum Monte Carlo approaches 
can be substantially damped by performing the sampling 
in each of the subspaces separately. Further, a rate de- 
scription developed for single particle tight binding sys- 
tems and known to be applicable whenever the dynamics 
is sufficiently incoherent, could be generalized to capture 
also the many body relaxation. This approximation is of 
great importance as it not only allows to obtain a deeper 
physical insight into the numerical QMC data, but also 
to investigate the long-time domain which is typically not 
accessible by the PIMC. 

We have illustrated these findings for one-dimensional 
Hubbard chains with sites of up to d — 4 and consider- 
ably long times, which may be also of experimental in- 
terest. In particular, small arrays of quantum dots would 
be appropriate candidates, where the above symmetries 
could directly be exploited to control dynamical features 
of two electron transport. In fact, a fully controllable de- 
vice consisting of two coupled quantum dots has already 
been realized 0]. An additional external magnetic field 
further opens a way to break or restore the fundamental 
spin-pcrmuation symmetry. Note that although here we 
have been mostly concerned with incoherent dynamics, 
the above findings also apply to the domain of coherent 
transport as shown for d = 3 in Ref. [Til ]. 

Systems with more than two electrons or systems with 
indistinguishable particles can be analysed along the 
same lines as above, however, with differences in detail. 
In this context the case of many boson systems coupled to 
a dissipative bath, i.e. a bosonic Hubbard model with dis- 
sipation, is certainly of great interest, e.g. for condensed 
ultra-cold atomic gases. Recently, we have already shown 
that for a two boson system a DFS exists only for an odd 
number of lattice sites (2(| . Work to extend this study to 
a full characterization of the symmetry properties as well 
as a generalization of the rate description are currently 
in progress. 
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APPENDIX A: DISCRETIZED PATH INTEGRAL 
FOR MANY-BODY SYSTEMS 

The derivation of a discretized version of the path- 
integral expression l|34|l is very similar to the correspond- 
ing case of a single particle (see, e.g., Refs plElHiEl)- 
Therefore, we here restrict ourselves mainly to details 
specific to the many-particle case. 

We begin with inserting 2q complete sets of projectors, 
one after each short-time step r = t/q, 

A{t) = tr{W(t)A} 

= E (IT /^) <si,^i|V^(0)|s 2 ,^ 2 ) 
I-Jl,|.jl<s \^= 2 / 

i=3 

x(s q ,x q \A\s q+1 ,x q+1 ) 

2q 

x II (sj-i,Xj-i\e 

j=q+2 



■irH/hl 



s*i. Xj 



(Al) 



where \x) is a suitable basis of the environment, and s = 
(s^, s^) denotes the system's coordinates with respect to 
^loc ©■ Exploiting the symmetric Suzuki splitting, 



= (s i _ 1 ,2f i _ 1 |e^+^)/ 2fi 



xe ^ S /V^+ ffB )/ 2 Vi,^) + 0(r 3 ) 
exp/^^Vi-T + ^-i-TlE 



x^lcxp - 



x exp< 



+ 0(r 3 ). 



IT 

2h 



-a(sj-i ■ e) 2J c a X a + H B 

a 



where 



(A2) 



(A3) 



denotes the free system's propagator over a time inter- 
val r, one obtains an expression where the environmental 
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degrees of freedom can be integrated out exactly due to 
the harmonic nature of Hb and Hj, leading to a dis- 
cretized version $[{sj}] of the influence functional (|35JI . 
Since, after replacing sj • e with the corresponding single- 
particle coordinate, the environmental terms in Eq. (|A2|) 
are then identical to those of the single-particle case (cf., 
e.g., Ref can be easily obtained from 

the single-particle expression without the need of redo- 
ing the tedious integrals. Rewriting the system's coor- 
dinates with the aid of the forward and backward paths 
introduced in Sec. IVI Al 



in 

3 

in 

3 



Ml 
j 

Ml 

b 2q+2-j 



for 1 < j < q + 1 , 

for q + 2 < j < 2q + 1 . 



(A4) 



the exponent of the influence functional becomes, in 
terms of the sum and difference coordinates ffj = (sj + 
SpS'j + sj) and ^ = (sj - sj, s] - sj), respectively, 

q 

3=2 



The real valued . . . Aj , Xj , and Xj are calculated accord- 
ing to 

A Q + iX Q = Q(t), 

Ki+iX l = Q((1-1)t)+Q((1 + 1)t)-2Q(It) 
for 1 < I < q - 2 , 
X, = Im{Q((j - 2)r) - Q((j - l)r)} 

for 2<j<q, (A6) 

where 

Q(t) = - [ du^^{coth(h/3uj/2)[l-cos(ujt)} 



+ i sin(cijt)} 



(A7) 



is the twice integrated bath-autocorrelation function 
L(t) 0- We thus arrive, as a discretized version of 
Eq. lEU), at 



( A (t)) = a ( f fl+ 1 ) Yl K ( f fj^jlVj+l,^3+^) 

+ E { ■ e j A i-* ' e j x exp{-*[{^,^.}]} , (A8) 

+ i (fj ■ ej Xj-k (ffk ■ ejj . (A5) which can now readily evaluated with PIMC schemes. 



T. Heinzel, Mesoscopic Electronics in Solid State Nanos- 
tructures ( Wiley- VCH, 2003). 

J. Gorman, D. Hasko, and D. Williams, Phys. Rev. Lett. [15 
95, 090502 (2005). 
J. Jortner and M. R. (eds.), Molecular Electronics (Black- [16 
well ScL, 1997). 

P. Hanggi and M. Ratner, Chem. Phys. 281 (2002). [17 
M. A. Nielsen and I. L. Huang, Quantum Computation [18 
and Quantum Information (Cambridge, 2000). 
A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. A. [19 
Fisher, A. Garg, and W. Zwerger, Rev. Mod. Phys. 59, 
1 (1987). [20 
U. Weiss, Quantum dissipative systems (World Scientific, 
1999), and references therein. [21 
A. Tsvelick and P. Wiegmann, Adv. Phys. 32, 453 (1983). [22 
A. Recati, P. O. Fedichev, W. Zwerger, J. von Delft, and [23 
P. Zoller, Phys. Rev. Lett. 94, 040404 (2005). [24 
J. Jortner and M. Bixon, eds., Electron Transfer-from 
isolated Molecules to Biomolecules, vol. 106/107 of Adv. [25 
Chem. Phys. (1999). 

L. Miihlbacher, J. Ankerhold, and A. Komnik, Phys. Rev. [26 
Lett. 95, 220404 (2005). 

F. Essler, H. Frahm, F. Gohmann, A. Klumper, and [27 
V. Korepin, The One-dimensional Hubbard Model (Cam- 
bridge University Pr., 2005). [28 
Y. Meir and N. Wingreen, Phys. Rev. Lett. 68, 2512 [29 
(1992). 

A. Ferretti, A. Calzolari, R. D. Felice, F. Manghi, M. J. 



Caldas, M. B. Nardelli, and E. Molinari, Phys. Rev. Lett. 
94, 116802 (2005). 

E. G. Petrov and P. Hanggi, Phys. Rev. Lett. 86, 2862 
(2001). 

F. Kaiser, M. Strass, S. Kohler, and P. Hanggi, Chem. 
Phys. 322, 193 (2006). 

R. Egger and C. Mak, Phys. Rev. B 50, 15210 (1994). 
L. Miihlbacher, J. Ankerhold, and C. Escher, 
J. Chem. Phys. 121, 12696 (2004). 

L. Miihlbacher and J. Ankerhold, J. Chem. Phys. 122, 
184715 (2005). 

M. Suzuki, Quantum Monte Carlo Methods in Condensed 

Matter Physics (World Scientific, 1993). 

J. Hirsch and D. Scalapino, Phys. Rev. B 29, 5554 (1984). 

J. Schulte and M. Bohm, Phys. Rev. B 53, 15385 (1996). 

J. Hirsch, Phys. Rev. B 65, 214510 (2002). 

S. Tornow, N.-H. Tong, and R. Bulla, cond-mat/0502276 

(2005). 

A. O. Caldeira and A. J. Leggett, Physica A 121, 587 
(1983). 

L. Miihlbacher, C. Escher, and J. Ankerhold, 
Chem. Phys. 322, 200 (2006). 

J. Sakurai, Modern Quantum Mechanics (Addison- 
Wesley, 1994). 

R. A. Marcus, J. Chem. Phys. 24, 966 (1956). 
S. Bose, Phys. Rev. Lett. 91, 207901 (2003). 



